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I. INTRODUCTION 


Since its introduction to the underwater acoustics 
community (Hardin and Tappert, 1973), the parabolic wave 
equation has stimulated a considerable amount of interest. 
The first solution programs used a split-step fast Fourier 
transform method to solve the parabolic equation; however, 
other solution techniques have been developed (McDaniel and 
Lee, 1982). One of the motives for developing alternative 
solution techniques is that problems arise when the Fourier 
transform encounters an interface between two media having 
diffarent sound speeds or densities (Lee and Botseas, 1982). 

One alternative solution technique uses an implicit 
finite-difference (IFD) solution method. The IFD method is 
unconditionally stable and has the capability to incorporate 
desired interface conditions. Implicit finite-difference 
methods for solving parabolic equations have been studied 
extensively by many authors. 

A computer program that utilizes the IF) method to solve 
the parabolic equation has been developed and is examined in 
detail in this thesis. The computer program predicts 
acoustic propagation loss in environments having user- 


specified bettom profiles. The program preserves continuity 


of pressure and continuity of the normal component of 


[a me - 
iv 7 **- 
= le AD 
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particle velocity at an interface between media having 
different sound speeds and densities. 

The program utilizes concepts developed by earlier 
authors. The use of the IFD method to solve the parabolic 
equation in underwater acoustics was developed by Lee and 
Papadakis (1979). The mathematical treatment of horizontal 
and sloping interfaces was developed by McDaniel and Lee 
(1982) and Lee and McDaniel (1983), And finally, the 
program utilizes some design features of an earlier computer 


program developed by Lee and Botseas (1982). 
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br II. PARABOLIC EQUATION 


A. INTRODUCTION 

The parabolic equation is an approximation to the 
: elliptical wave equation. The derivation of the parabolic 
: equation begins with the reduced wave equation (Helmholtz 


equation) in the form 


7 vy? p+ k? = 0 2.1 

4 

“ or 

“ 

it 

vy? p+ ko2n2p = 0 eure 
7 where 
3 k # wave number (= w/c) 

: k, = reference wave number (= w/c.) 


n = index of refraction (= c,/c) 
i S p # time independent factor of complex pressure 
c »# sound speed 
Cg = reference sound speed 


w = angular source frequency (= 27 f) 


id For the case of cylindrical symmetry (2.2) becomes 

§ Pry + (1/r) Dy, + Deg + k,?n@p = 0 2.3 
a 

It is then assumed that p is of the form 

4 


p = u(r,z) S(r) 


- os wf 


where u is a function of both range and depth and S is a 


P.. function of range only. Substitution of (2.4) into (2.3) 
a 

¥ and separation of variables shows that S(r) must satisfy 
x 

a: 

“4 11 
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=> 
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can 
fi Bessel's equation of zero-order. For exp (-iwt) time 
“ 

oA 


dependence and outgoing waves, the solution is the zeroth- 


» ead order Hankel function of the first kind, 
i s(x) = HA() (Kr), 
a ° fe) 

, 

: Further, u(r,z) must satisfy 


_— 1 2 
i Upy + Uge + (- * = SZ) ur + ko(n? - 1) uz 0. 2.5 
r Ss 


With the help of the far-field asymptotic approximation for 


FE the Hankel function, (2.5) can be reduced to 

: eS 

FF Ure + Uz, + 2ikguy + ko2(n* - 1) u = 0. 2.6 
PN 
Ln We now assume that u varies slowly with respect to range, 


; 
} 
F: |e | &< 2 kot | 27 
combining (2.6) and (2.7) results in 

| 


= 

7 Up, + 2ikou, + ko2(n2 - 1) u = 0. 2.8 
& Rearranging (2.8) results in the parabolic equation in the 
Me form 
[ nes up = a(ky,r,z) u + b(k,,r,z) uz, 2.9 
x where 2 
tr a(ko,r,z) = (ik,/2) ([n“(r,z) - 1] 
oN 
Le b(ko,r,z) = i/2k,. 
oe The assumption (2.7), fundamental to the parabolic equation 
es 

a method, is equivalent to neglecting back-scatter. 


B. SPLIT-STEP FAST FOURIER TRANSFOKM SOLUTION 
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[8 For the first few years after its introduction into 
a the acoustical community the parabolic equation was sclved 
a a, 
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exclusively with the help of the split-step fast Fourier 
transform (SSFFT) method developed by Tappert and Hardin 


(Jensen and Krol, 1975). In this method, u in (2.8) is 


Zz 
represented by the inverse transform of its Fourier 
transform. The SSFFT method requires periodic boundary 
conditions in z because of the finite Fourier transform. 
This is handled by introducing an artificial, horizontal, 
pressure release bottom below the physical bottom. The 
SSFFT method is unconditionally stabie (Brock, 1978). 

The SSFFT method has been implemented by Jensen and 
Krol and by Brock. Detailed descriptions can be found in 
publications of Jensen and Krol (1975) and Brock (1978). 

2. Interface Treatment 

Errors introduced by the SSFFT method are 
proportional to the range step and to No, where n is the 
index of refraction (Jensen and Krol, 1975). Because the 
index of refraction has a large change across the water- 
sediment interface n,, will be large and thus the error will 
be large. To reduce this error, a very small horizontal 
range step must be used. However, this results in very long 
computer execution time. The problem of a discontinuity in 
sound speed at the water-sediment interface and the 
resultant difficulties in solving shallow water propagation 
problems using the SSFFT method are addressed in Jensen and 


Krol (1975). 
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Another more serious problem with the SSFFT method 


is that it neglects any density difference between the water 


‘Oy and the sediment. A density difference can be important in 
ea . that it influences the reflection coefficient. The problem 
becomes more significant as the density difference becomes 
* larger. 

% In summary, the discontinuities in sound speed and 
ot in density at the water-sediment interface cause problems 
ni for the SSFFT method. The SSFFT method is therefore 
a intrinsically better suited for deep water propagation 
= environments for which the water-sediment interface is an 
; unimportant feature. 

q 

. C. IMPLICIT FINITE-DIFFERENCE SOLUTION METHOD 

oH .In 1979 Lee and Papadakis introduced the Crank-Nicolson 
ri ; implicit finite-difference method to solve the parabolic 

a equation for underwater acoustic propagation. The Ccank- 
a Nicolson method uses a second-order central difference 
% formula to approximate u,, in (2.9) and casts the problem in 
Er the form of a tridiagonal matrix system. A representative 
% row in the matrix system (the m™" row) is 
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i k = horizontal range increment 
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h = vertical depth increment 
Supercripts are used to indicate range indices and 
subscripts are used to indicate depth indices. In (2.10) 
the field is known at range index n and is to be solved at 
range index n+ 1. Therefore, the right hand side of (2.10) 
reduces to a single, known value and the solution field is 
advanced from range index n to range index n + 1 by solving 
the tridiagonal system of equations. 

The IFD scheme is consistant, unconditionally stable and 
it converges to the theoretical solution as the range and 
depth increments tend to zero (Lee et al., 1981). An 
advantage of selecting an implicit scheme over an explicit 
scheme is that an explicit scheme is only conditionally 
stable (Lee and Papadakis, 1979). Another advantage of the 


implicit scheme is smaller errors. More detailed 


Seti ts gM A Ot ge rN ea iy Ses sae ee 


8 information addressing both implicit and explicit solutions 
ie of parabolic equations can be found in Gerald (1980). 

a The first IFD scheme handled discontinuities in the 
« speed of sound profile but did not consider the effects of 
Re r density discontinuities. It therefore did not correctly 
gy treat the interface between media having different 
oe densities. 

" D. IMPLICIT FINITE-DIFFERENCE METHOD: TREATMENT OF A 
& HORIZONTAL INTERFACE 

“a In 1982 McDaniel ard Lee introduced a scheme for 
x incorporating a horizontal interface into the IFD method. 
x The interface separates two media with different sound 
BN speeds and densities (Figure 1). 
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- Figure 1. IFD Treatment of a Horizontal Interface 


ty 
Ss 
ra 16 


The scheme preserves continuity of pressure and continuity 


i of the normal component of particle velocity across the 
H interface and does not affect the stability of the IFD 
method. 


The interface equation in the tridiagonal matrix system 
that results from incorporating continuity of pressure and 
continuity of the normal component of particle velocity is 


(McDaniel and Lee, 1982) 
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f°; = density in layer 1 (water) 
/°2 = density in layer 2 (sediment) 
a and b are defined in (2.9) 
Incorporating the horizontal interface into the IFD method 
requires inserting (2.10) for the row in the tridiagonal 


matrix system that corresponds to the interface. 
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The error in the solution is 
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eerie y. 


O (k? + kh) 
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on the interface and 

© (k? + kh?) 
in a continuous medium (McDaniel and Lee, 1982). 

E. IMPLICIT FINITE-DIFFERENCE METHOD: TREATMENT OF A 

SLOPING INTERFACE 

In 1983, Lee and McDaniel extended their treatment of an 
interface between two media to include the case of a sloping 
interface. As for the case of a horizontal interface, the 
treatment of a sloping interface preserves continuity of 
pressure and continuity of the normal component of particle 
velocity at the interface between media having different 
densities and sound speeds. 

The problem of a sloping interface is separated into two 
cases: downslope, and upslope. Each case requires inserting 
two new rows into the IFD tridiagonal matrix system. One 
new row is required at the level corresponding to the 
interface at the range where the solution is known and the 


second new tow is required at the level corresponding to the 
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interface at the range at which the solution is to be 
solved. Therefore, four new equations are required to 

cover both the downslope and upslope cases. The four 
sloping interface equations are derived and shown in Lee and 
McDaniel (1983). The equations are somewhat involved but 
they are of the same tridiagonal form as the original IFD 
matrix equations. The error for points on or adjacent to a 
sloping interface is (Lee and McDaniel, 1983) 

O (k? + kh). 
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III. COMPUTER IMPLEMENTATION 


A. INTRODUCTION 

An IFD solution program that implements (2.10), (2.11) 
and the four sloping interface equations has been developed 
to solve underwater propagation problems. The program is 
written in FORTRAN using single precision, complex 
arithmetic and has been installed on the Naval Postgraduate 
School's IBM-3033 digital computer. Appendix A contains a 
program listing. 

The solution program consists of one main program and 20 
subroutines. A modular program construction was selected 
for flexibility and clarity. 

Within each routine a structured programming approach is 
utilized. The structured program format, coupled with 
generous commenting, makes the program relatively easy to 
trace through. 

As presently installed on the NPS computer the solution 


program is run interactively. Appendix B contains details, 


B. GENERAL CHARACTERISTICS 

The solution program handles the following environmental 
conditions: range independent sound speed profile in the 
water column, range dependent bottom profile and iso-speed, 


iso-density sedimentary bottom layer. The program utilizes 
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a Gaussian starting field and an artificial pressure release 
z:. surface in the sediment at a user-specified depth. (An 
d artificial pressure release surface is not required for the 


IFD solution method; however, such a surface permits 


an efficient solution for the pressure field when it is known 
Ry to tend to zero at great depths in the botton. ) 
Ee Attenuation in the water and sediment is introduced 
a using complax indices of refraction. Artificially strong 
ine attenuation is applied in the lower portion of the sediment 
BS layer to enhance attenuation of the field above the 
a artificial pressure release surface. 
me A user-specified bottom profile is input as a series of 
Ss linear bottom segments. The range step along a horizontal 
sts bottom segment is set to a user-specified value. The range 
ao step along a sloping bottom segment is automatically set by 
Ps ‘ the program so that the sloping bottom intersects the next 
ae vertical grid point. Bottom modifications are required in 
® certain situations to meet the requirements that the range 
Ss step not be too large and that the interface must pass 
ah through a grid point at every range at which the pressure 
ee field is solved. For a very gently sloping bottom, if the 
of calculated range step exceeds a user-specified value then 
= the program will automatically model the bottom as a series 
rhe of level and sloping sections. For these cases, the 
aa difference between the modified bottom and the user- 
ba specified bottom is always less than or equal to one-half 
Fes 
A 
tet 21 


the vertical grid increment. The user is informed if the 
bottom is modified. 

It is foreseen that future enhancements will increase 
the program's generality. In particular, changes to allow 
range dependent sound speed profiles, sound speed profiles 
in the sediment and a user-specified starting field should 
be relatively simple. The modular construction of the 


program facilitates these types of changes. 


C. MAIN PROGRAM IFD 

IFD is the main program. It controls program execution 
and calls subroutines as appropriate. 

The first executable statement in IFD calls subroutine 
ERRSET, a system subroutine peculiar to the NPS computer 
that correctly sets a variable value to zero when an 
underflow condition exists. Most computer systems do this 
automatically; however, depending on the particular system a 
call similar to ERRSET may be required. 

IFD then calls subroutine READ to read input data, 
subroutine SVPW to calculate the sound speed eat grid points 
in the water column, subroutine INITAL to initialize 
constants and variables, and then subroutine MATCON to 
calculate matrix constants. Subroutine SFIELD is then 
called to calculate the Gaussian starting field, followed by 
subroutines WRITE1 and PRINT! which write and print output 


data. In the context of this program, "write" refers to 
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writing unformatted data into a file to be used by the 
plotting routine and “print” refers to writing formatted 
data into a file which can be sent directly to the printer. 
IFD then calls subroutine NEWSEG which is the beginning 
of a loop that is called every time a new linear bottom 
segment is reached. NEWSEG calculates variables that 
characterize a new bottom segment. The next call is to 
subroutine NEWMAT which calculates matrix elements for the 
new battom segment and advances the solution field one range 
step. IFD then enters a loop that advances the solution one 
range step for every pass through the loop. Inside the loop 
the range markers are advanced and the solution is advanced 
one step for the downslope, level, upslope, modified bottom 
downslope or modified bottom upslope situation as 
appropriate. In addition, the artificial attenuation 
mentioned earlier is applied by calling ATTENU and calls are 
made to WRITE2 or PRINT2 as required. Finally the range is 
checked to see if it has advanced to the maximum range 
specified tor the problen. If it has, then IFD calls 
subroutine END which passes appropriate messages to the 


terminal and stops program execution. 


D. SUBROUTINES 
1. Subroutine READ 
Subroutine READ is called by IFD to read input data 


from unit number NIU = 51. Input data are read in free 
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format and data are transferred back to main program IFD via 
common blocks. REAC contains error checks for (1) input 
data insufficient and (2) the final depth in the sound speed 
profile unequal to the maximum depth in the water column. 
If either of these error conditions exists, READ issues an 
appropriate error message to the terminal and stops 
execution, 
2. Subroutine SVPW 

Subroutine SVPW calculates the vertical grid spacing 
used in the water and sediment. It also calculates the 
speed of sound at each of the grid points in the water 
column. Linear interpolation is used to calculate the sound 
speed at grid points between points on the user-specified 
sound speed profile. 


3. Subroutine INITAL 


Subroutine INITAL initializes constants and 
variables. If the user inputs 0.0 for the value of the 
reference sound speed then {NITAL sets the reference sound 


speed c. to the sound speed averaged over the deepest water 


fe) 
column. If the user inputs 0.0 for the value of the maximum 
range step then INITAL sets the maximum range step to the 
reference wavelength, 

DRMAX = XLAMDA . 
Setting the maximum range step to the reference wavelength 


is somewhat arbitrary; however, until the actual limit on 


the range step is better understood it serves as a rough 
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rule of thumb. Finally, if the user inputs 0.0 for the 
value of the range step along a horizontal interface then 
INITAL sets the range step to half of the reference 
wavelength, 
DRLVL = 0.5 * XLAMDA. 
The default range step along a horizontal interface is haif 
the default maximum range step. 
4. Subroutine MATCON 

Subroutine MATCON calculates constants needed to 
compute tridiagonal matrix elements. Most of the constants 
computed in MATCON have no direct physical significance but 
contribute to computational efficiency. Attenuation in both 
the water and sediment is calculated with the help of a 


complex index of refraction n, 


Cy BETA 
n={— (1 + L —— ea } 
Cy 54.575054 
or 
2 2 
Cc c BETA 
n2= acd + L Zn a 
C4 Cy 27 .287527 
where 


BETA = attenuation (dB/wavelength) 
Co = reference sound speed (m/s) 
Cy = sound speed (m/s) at point j 
54.575054 = conversion factor used in converting 


db/wavelength to nepers/meter. 
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5. Subroutine SFIELD 
Subroutine SFIELD calculates the Gaussian starting 
field at range r = 0. This subroutine is identical with 
that of Lee and Botseas (1982); both yield the starting 
field suggested by Brock (1978), 
U(I) = CMPLX (PR, 0.0) 


zM-zs | 2 -z7M-zs |2 
GW GW 


PR = GA [e -e ] 


where 


RQ 
K 
a 


depth (m) of grid point 
2S = source depth (m) 
GW 


Gaussian width (m) (2 2/FK) 
FK 


reference wave number (1/m) (# 2 f/c,) 


GA = Gaussian amplitude [= (2/FK) 1/2 /qw 


6. Subroutine WRITE! 

Subroutine WRITE1 outputs data to a file that is used 
by the plotting routine, This output file corresponds to 
unit file number NOU = 52. 

7. Subroutine PRINT! 

Subroutine PRINT1 outputs formatted data toa file 

that can be sent to the printer. The output file 


corresponds to unit file number NPOUT = 55. 
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8. Subroutine NEWSEG 
i Subroutine NEWSEG is called at the start of each new 
linear bottom segment. NEWSEG computes and initializes 
. variables that depend on characteristics of the segment. 
Bal One of the variables initialized is ISLOPE which is a slope 
flag having value 1 if the bottom slopes down, 2 if the 
bottom is horizontal, 3 if the bottom slopes up, 4 if the 
bottom slopes down but must be modified because the slope is 
very small, or 5 if the bottom slopes up but must be 
& modified because the slope is very small. NEWSEG also 
2 


issues error or warning messages as appropriate. 


9. Subroutine NEWMAT 

Subroutine NEWMAT calculates matrix elements for the 
X and Y matrices. The Y matrix corresponds to the range at 
i which the solution field is known and the X matrix 
: corresponds to the range at which the solution field is to 
be found. The Y matrix is multiplied by the known solution 
a! field to obtain the right-hand side column vector needed to 
® solve the tridiagonal system. 

NEWMAT sets up the tridiagonal matrix system for the 
new bottom segment and than calls TRIDG to solve the system 
at the first range step. It then calls ATTENU to apply 


artificial attenuation, calls WRITE2 or PRINT2 as required, 


cit Pe 


and finally updates the interface pointer that indicates the 


index of the grid point at the water-sediment interface, 
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10. Subroutine WRITE2 
Subroutine WRITE2 is basically a continuation of 
subroutine WRITE1. It outputs data to a file corresponding 
to unit file number NOU = 52. This is the file that is 
used by the plotting routine. At range intervals specified 
by the user, WRITE2 outputs range, depth, and the value of 
u(r,z). 
11. Subroutine PRINT2 
Subroutine PRINT2 is basically a continuation of 
subroutine PRINT1. It outputs formatted data to a file 
corresponding to unit file number NPOUT = 55. This file can 
be sent directly to the printer. At range and depth 
intervals specified by the user PRINT2 outputs tabular 
values of transmission loss and u(r,z). 
12. Subroutine TRIDG 
Subroutine TRIDG solves a linear tridiagonal matrix 
system. TRIDG is a modified version of subroutine TRID as 
list din Gerald (1980). The major modifications to 
subroutine TRID involved introducing arrays CTWO and CR to 
preserve the original matrix element values and to make the 
routine more efficient. Introducing the two new arrays 
requires more storage space but results in a substantial 
savings in execution time, particularly for the case cof a 


horizontal interface. 
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13, Subroutine TRIDL 
Subroutine TRIDL is a modified version of subroutine 
TRIDG. TRIDL differs from TRIDG in that it does not compute 
arrays CTWO and CR but rather uses the array values 
calculated in TRIDG. 
14. Subroutine DOWN 
Subroutine DOWN updates the tridiagonal matrix and 
calls subroutine RHS to update the right-hand side for the 
case of a downward sloping interface. DOWN then calls 
subroutine TRIDG to solve the tridiagonal system of 
equations and finally updates the interface pointer. 
15. Subroutine UP 
Subroutine UP performs exactly the same tasks as 
subroutine DOWN, but for the case of an upward sloping 
interface. Subroutine TRIDG is again called to solve the 
system. 
16, Subroutine LEVEL 
Subroutine LEVEL is called to advance the solution ‘ 
for the case of a horizontal interface. For this case the 
tridiagonal matrix elements at the advanced range need not 
be changed from the previous calculation. Therefore, LEVEL 
need only update the right-hand side by calling RHS and 
then solve the system by calling TRIDL. 
17. Subroutine RHS 
Subroutine RHS computes the right-hand side of the 


tridiagonal system by multiplying tridiagonal matrix Y by 
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bE | 
a the known solution field U(I). The resultant right-hand 
= side column vector is stored in C(I,4). 

- 18. Subroutine SSLOPE 

BS Subroutine SSLOPE is call:d to advance the solution 
= in the case where the bottom has been modified. SSLOPE 
gan determines which of three cases a particular section falls 
Py into: a level section following a level section, a level 
me section following a sloping section, or a sloping section. 
™ For the case of a level section following a level section 
aN SSLOPE calls LEVEL to advance the solution. For the case of 
oe a level section following a sloping section SSLOPE updates 
| appropriate matrix elements, calls RHS and then calls TRIDG 
i to advance the solution. And in the case of a sloping 
ei ’ section SSLOPE updates matrix elements and calls either DOWN 
“a or UP as appropriate. 

ea 19, Subroutine ATTENU 

ey Subroutine ATTENU applies artificial attenuation to 
“4 the bottom portion of the sediment layer as suggested by 
& Brock (1978). The artificial attenuation matrix ATT(1) is 
ea calculated in subroutine NEWMAT. 

io 20. Subroutine END 

Be Subroutine END is called when the solution field has 

a reached the maximum range specified. END sends appropriate 

a messages to the terminal and stops execution. The messages 
3 are applicable to the program as installed on the NPS 
ney 
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computer but may not be appropriate for the program if 


we installed on another systen. 


E. INPUT DATA 


1, Input File 

“ The input data must be stored ina file 
¥ corresponding to unit number NIU as assigned in subroutine 
a READ. In its present form READ sets the input unit number 

to NIU = 51. If the user prefers to read the data from a 

¥ different unit (for example, a card reader), then variable 
x NIU in READ should be set equal to the appropriate unit 
, : number. 

2. Input Format 

5! The input data is read in free format. The input 
| card images (or input cards) are arranged as follows: 

A CARD CONTENTS 

¥ 1 FRQ, 2S, ZR, CO, N 

; where 

x FRQ = frequency (Hz) 

r ZS = source depth (m) 

a ZR = receiver depth (m) 

F (program will reset to depth of nearest 
_ grid point) 

CO # reference sound speed (m/s) 
: s | If CO = 0.0, CO is set to the sound speed 
: 4 ’ averaged over the daepest water column. 
FF 
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N = number of vertical grid points 


R If the user desires that every integer 
~ depth value correspond to a grid point then 
BS f (neglecting dimensions) N should be set to 
FS an integer multiple of ZLYR2, the depth of 
the pressure release surface. 
% 

= CARD CONTENTS 

4 2 RMAX, DRLVL, DRMAX, WDR, PDR, PDZ 

s where 

: RMAX = maximum range (m) of solution 

: DRLVL = range step (m) for marching solution 
o along horizontal interface 

3 If DRLVL = 0.0, then DRLVL is set to 1/2 


wavelength. 
If DRLVL is greater than DRMAX, then 


DRLVL is set to DRMAX, 


iy ede EalS 


DRMAX = maximum allowable range step (m) 


If DRMAX = 0.0, then DRMAX is set to 1 


Ys wavelength. 
° WDR =z range increment (m) at which solution is 
4 written to file used by plotting routine 
: PDR z range increment (m) at which solution is 
- printed 
i PDZ = depth increment (m) rounded to nearest DZ 
> at which solution is printe 
. t which solution is printed 
: 32 
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a CARD CONTENTS 


S 3 BR(1) , BZ(1) 

if 4 BR(2) , B2(2) BOTTOM PROFILE 

3 ; . S§5 BR(3) , B2(3) Range and depth of water (m). 
8 ; . . Maximum number of points = 100, 
4 ? 7 P Program will reset depths to 

: . : ° nearest grid point. 

zi N , 

z N+1 -1 -1 This card marks end of bottom 
y profile. 

x N+2 ZLYR1, RHO1, BETA! 

7 where 

x ZLYR1 «= maximum water depth (m) 

é RHO1 = density of water (g/cm?) 

BETA1 = attenuation in water (dB/meter) 

can If BETA1 is less than 0.0 then program 
. calculates BETA1 with an empirical 
r formula (Brock, 1978). 

a N+3 ZSVP(1), CSVP(1)| SOUND SPEED PROFILE 

N+4 ZSVP(2), CSVP(2) Depth (m}) and sound speed 
* : . ‘ (m/s). 

RS : . ZSVP(1) must equal 0. 

3 N+M The last depth must equal 
Sa ZLYR1. 
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CARD CONTENTS 


N+M+1 ZLYR2, RHO2, BETA2, C2 


where 
ZLYR2 = 
RHO2 2 
BETA2 = 
C2 = 


N+M+2 ZABLYR 


where 


ZABLYR = 


F. PROGRAM OUTPUT 


depth (m) of pressure release surface at 
bottom of sediment layer 

density (g/cm?) of sediment 

attenuation (dB/wavelength) in sediment 


sound speed (m/s) in sediment 


depth (m) of upper surface of artificial 
attenuation layer in sediment. 


ZABLYR should be about 3/4 of ZLYR2. 


1. Output Printer File 


The program outputs formatted data to a file 


corresponding to unit number NPOUT which is set to 55. This 


formatted data file may be sent to the printer if desired. 


On another system the user may elect to assign NFOUT to the 


unit number corresponding to the printer and thereby send 


the formatted data directly to the printer. 


2. Qutput Plotter File 


The program outputs unformatted data to a file that 


is used by the plotting routine. The unit number for this 


file is NOU which is set to 52. The output data in this 
file are stored as follows: 


LINE CONTENTS 


RMAX = maximum range (m) of solution 


2 RA, 2R, U 
3 RA, 2R, U 
where 
RA = range (m) 
ZR = depth (m) 


U = complex u at specified range and depth 
3. Terminal Output 
Certain WRITE statements in the program specify unit 
number 6. Unit number 6 on the NPS computer for an 
interactive program corresponds to terminal output. If the 
program is not run interactively then WRITE statements with 
unit number 6 may be deleted. Any pertinent information 
passed to the terminal during an interactive run is also 


passed to unit number 55, 


G. PLOTTING PROGRAM 
Appendix C contains a listing of the IFD plotting 


program installed on the NPS computer. The filename and 
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filetype of the program are PLOTIFD FORTRAN. This program 
was written separately from the IFD program because it was 
recognized that different computer installations have 
different plotting facilities. For users with different 
facilities the program will be a helpful raference. 


The program reads data from unit number NOU = 52 which 


corresponds to a file with filename/filetype IFDOUT PLOTTER. 


Details concerning using the plotting routine are included 


in Appendix B. 


IV. TEST CASES E 


A. HORIZONTAL INTERFACE CASES 
The IFD program treats a horizontal interface using the 

same theoretical approach as the IFD program published by 
| Lee and Botseas (1982). Throughout the remainder of this 
‘ thesis the Lee and Botseas (1982) program will be called the 
: FINITE-DIFFERENCE program and the program presented in this 
thesis will be called the IFD program. Two cases were run 
fa to confirm that the IFD program is in agreement with the 
FINITE-DIFFERENCE program for horizontal interfaces. 


¥ 1. Isospaed Shallow Water 


This case, first published by Jensen and Kuperman a 
(1979), considers propagation in a shallow water, isospeed 
environment. The water depth is 100 meters and the solution 
field is calculated out to 25 kilometers. The sound speed 
f: is 1500 m/s in the water and 1550 m/s in the sediment. 
Density and attenuation in the sediment are 1.2 g/cm? and 1 
N dB/wavelength respectively. The source and receiver are 
both at 50 m and the source frequency is 500 Hz. 

Solutions obtained using a normal mode program 
(SNAP), a split-step fast Fourier transform program (PAREQ) 
and the FINITE-DIFFERENCE program are shown in Figure 2. 
SNAP and PAREQ are programs that were developed at SACLANT 


Centre and are discussed further in Jensen and Kuperman 
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(1979). The solution obtained using the IFD program is 


shown in Figure 3. All of the solutions are in excellent 
agreement. 
The input runstream that produced the results shown 


in Figure 3 for the IFD program is as follows: 


Input Runstream 


500 50 50 0 500 

25000 5 5 50 5000 50 
0 100 

25000 100 

-1 -1 

100 1.0 -1.0 

0 1500 

100 1500 

250 1.2 1.0 1550 

200 


2. Horizontal Interface 

This case, called the “horizontal interface problem" 
in Lee and Botseas (1982), considers propagation in an 
environment with the sound speed profile shown in Figure 4, 
Scurce frequency is 100 Hz, source depth is 30 m and receiv- 
er depth is 90 m. The density in the bottom is 7.1 g/cm? 
and the sound speed in the bottom is 1505 m/s. No attenua- 
tion is applied in the water or sediment using complex 
indices of refraction; however, artificial attenuation is 


applied in the lower portion of the sediment layer. 
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‘e Figure 4. Horizontal Interface Case 


The solution obtained using the IFD program is shown 
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in Figure 5. This solution is in excellent agreement with 


the FINITE-DIFFERENCE solution shown in Lee and Botseas 


PL stir 


(1982). 
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The input runstream for this case is as follows: 
Input Runstream 
100 30 90 0 600 


0 240 
20000 =. 240 

A “| 

240 1.0 0.0 
0 1500 

120 1498 

240 1500 


1200 21 0.0 1505 
512 


B. RANGE-DEPENDENT CASES 

The following range-dependent cases were solved by 
Jensen and Kuperman using SNAP, the normal mode program, and 
PAREQ, the SSFFT program (Jensen and Kuperman, 1979). The 
cases were also solved by Lee and Botseas using the FINITE- 
DIFFERENCE program (Lee and Botseas, 1982). The FINITE- 
DIFFERENCE program treats the sloping interface as a "stair 
step" and uses the interface conditions appropriate for a 
horizontal interface. The IFD program handles the sloping 
interface using the interface treatment derived by Lee and 


McDaniel (1983). 
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es 
i 1. Deep-to-Shallow Water 
- ides Heer tosses. 

ee This case considers propagation in an environment as 
is depicted in Figure 6. The problem is solved for a bottom 
ae with a 8.5 degree upslope and one with a 0.85 degree 
4 

34 

oN upslope. Source frequency is 25 Hz, source depth and 
re receiver depth are 25 meters. Sound speed in the water is 
“nn 

os 1500 m/s. In the sediment, sound speed is 1600 m/s, density 
ve 

bE . is 1.5 g/cm? and attenuation is 0.2 dB/wavelength. 
s 0 0 20 30 40 km 

n ~ 80 
(Wal z 100 C,* 1300 m/s 

i = 150 Cy * 1600 m/s 

me a, 200 Peis Yom? 

FS Ww 250 pg"! 

pase! © 300 Be° 0.2 dB/Xd 

ie 350 

. Figure 6. Deep-to-Shallow Water Case 

" 

ay 

o The results for the 8.5 degree upslope case as 
- produced by SNAP,PAREQ, and FINITE-DIFFERENCE are shown in 
‘my Figure 7. The results as produced by IFD are shown in 

4 
; £4 Figure 8 The difference between the results produced by 
et SNAP and PAREQ is attributed to failure of the "adiabatic" 
4, 

int 

ey theory underlying the SNAP program (Jensen and Kuperman, 
a 

Pe) 

a, 1979). As determined from the input runstream the FINITE- 
Ka! DIFFERENCE results were obtained using 1.0 rather than 1.5 
ae 

ae g/cm3 for the density in the sediment (Lee and Botseas, 
4" 

yf 

ae 1982). 
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The input runstream for the IFD program that 
produced the results shown in Figure 8 is as follows: 


Input Runstream 


25 25 25 0 1000 
40000 10 ) 100 10000 25 
0 350 

10000 350 

12000 50 

40000 50 

= a 

350 1.0 = 

0 1500 

350 1500 

1000 1.5 0.2 1600 

750 


Appendix D contains the printed output produced by the IFD 
program using the above runstream. 

The results for the 0.85 degree upslope case as 
produced by SNAP and PAREQ are shown in Figure 9. The 
results produced by IFD are shown in Figure 10. 

2. Shallow-to-Deep Water 

This case considers propagation in the environment 
depicted in Figure 11. The environment is exactly the same 
as the deep-to-shallow water environment except that the 
shallow and deep portions have been reversed and thus the 


bottom slopes down rather than up. 
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Figure 9. Propagation Loss Versus Range for Deep-to-Shallow 
c, Water Case, 0.85 Degree Slope; SNAP and PAREQ 
Results 

The results for the 8.5 degrees downslope case as 
8 produced by SNAP and PAREQ are shown in Figure 12. As 
| ‘ before, the difference between SNAP and PAREQ is attributed 
x to failure of the SNAP program. The results produced by IFD 
i. are shown in Figure 13. 

The results for the 0.85 degree downslope case are 


shown in Figures 14 and 15. 


3. Comments 

& Differences between the results obtained using the 
* SNAP and PAREQ programs for the range-dependent cases are 
FS : discussed in Jensen and Kuperman (1979). The major 
‘ differences are attributed to the violation of the 
+] adiabatic assumption in the SNAP program. 
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Figure 11. Shallow-to-Deep Water Case 


boi m , 


60 --7-"- SNAP 


cmt 


PROPAGATION LOSS (d@) 


eo 


fe 


=F 


. RANGE (km) 


ees 


Figure 12. Propagation Loss Versus Range for Shallow-to- 
Deep Water Case. 8.5 degree Slope; SNAP and 
=” PAREQ Results 
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Figure 14. Propagation Loss Versus Range for Shallow-to- 
Deep Water Case, 0.85 degree Slope; SNAP and 

PAREQ Results 
Figure 16 shows the results produced by IFD for the 
8.5 degree, deep-to-shallow water case if 1.0 rather than 
1.5 g/om? is used for the density of the sediment. The IFD 
results obtained using 1.0 g/em? are in very close agreement 
with the PAREQ results shown in Figure 7. However, when the 
correct value of 1.5 g/cm3 is used the slope of the 
propagation loss curve beyond 15 km is less steep (Figure 8) 
and the results do not agree as well with the results 
produced by PAREQ. The observation that the slope of the 
propagation loss curve becomes less steep when 1.5 g/cm? is 
uged is qualitatively consistant because a higher density 


difference means a higher reflection coefficient which in 


52 


s}Insey aggI fedorts 
serbep sg°o ‘asep rz0eqemM deoq-03-MoTTeUS 
ioj obuey snszeq ssoT uotqzebedorzqg “Ss, eanbty 


(WA) SONUY 


0°OF 0°O£ 0°02 0°OT 0°0 

oa 

[ a} 

oO 
oe 

= 6 

aan |) 

o7 
6) 
a) 
a, a 
—" 
oO 

oz 

ae Sa 
oO 
T)) 
Ww 

~) — 

[o0 

o WwW 

yi 

= 

So 

. . a 
SA! EOS. CLEA? Beas | MPUTCe Sore 


HEN ERLE SSE Geet RAMAAARL FSO LES COREE 


A 

- turn results in more energy confined within the water layer. 

§ The difference in slope between the PAREQ and IFD results is 
attributed to the PAREQ program's apparent failure to 
¥ 
§ _ account for the affects of the density discontinuity at the 

# water-sediment interface. Because the IFD program correctly 

‘ a 

" accounts for density discontinuities the results produced by 
7 IFD are believed to be more accurate than those of PAREQ 
“4 

7 when interface interaction is important. 
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V. COMMENTS AND CONCLUSIONS 


The IFD method is an efficient, stable method for 
solving the parabolic equation. Use of the IFD method is 
particularly advantageous in shallow water environments 
where the water-sediment interface is an important 
parameter. 

The IFD program presented in this thesis incorporates 
continuity of pressure and continuity of the normal 
component of particle velocity across horizontal and sloping 
interfaces. The program's capability to incorporate the 
exact interface conditions on a sloping interface, to 
automatically determine step-size, and to modify the bottom 
as required for the case of a very gently sloping bottom are 
important features. 

Projected program enhancements include wide angle 
propagation (Lee and Gilbert, 1982), range-dependent sound 
speed profiles in the water, range-dependent sound speed 
profiles in the sediment layer, and multiple sediment layers 
with horizontal or sloping interfaces. These enhancements 
are listed in their approximate order of importance. The 
program's modular construction and structured style will 


facilitate implementation of these enhancements. 
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APPENDIX A: IMPLICIT FINITE-DIFFERENCE PROGRAM LISTING 
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APPENDIX B 


RUNNING THE IMPLICIT FINITE-DIFFERENCE 
PROGRAM ON THE NPS COMPUTER 
A. INTRODUCTION 
This appendix describes one simple procedure for running 


the IFD program on the NPS computer. 


B. COPYING FILES ONTO USER'S DISK 
Four files are needed to run the IFD program; the 


filenames and filetypes are 


IFD FORTRAN 
IFD EXEC 
PLOTIFD FORTRAN 
PLOTIFD EXEC 


The files are available from a computer account maintained 
by the underwater acoustics curriculum. To link with this 
account and obtain copies of the files the user should 
proceed as follows: 

(1) Log on terminal 

(2) Enter: CP LINK 0160P 191 195 RR 

(3) When prompted for read password enter: UX 

(4) Enters ACC 195 C 

(5) Enter: COPYIFD 
At this point the four files should reside on the user's A 
disk. 
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C. RUNNING THE IFD PROGRAM 

Before running the IFD program the user should define 
additional storage space, compile the program, and set up 
the input data file. To define additional storage space, 


(1) Enter: DEF STOR 1M 
(2) Enters: I CMS 


These two commands need only be entered one time for each 
terminal session; the storage space will remain for the 
entire session. 

To compile the program, 

Enters: FORTGI IFD 
The program need be compiled only one time unless the 
program is changed, in which case the new version should be 
recompiled. 

The final step before running the program is setting up 
the input data file which has filename and filetype IFDIN 
DATAIN. The user must create or modify this file so that it 
contains input data as described in Section III.E.2 of this 
thesis. For more information concerning how to create or 
modify files see NPS Technical Note TN-VM-02 which is 
available in the computer consultant's office. 

If the above steps are accomplished, the user can then 
run the program with 

Enters: IFD 
Shortly after entering this command the user will be 


prompted for a run identification. The run identification 
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ds an arbitrary identification label that will appear on the 
output printer file. Enter run identification as desired. 
At the end of the run the user is informed of the two 
output data files generated by the program. If desired by 
the user the output printer file (IFDOUT PRINTER) may be 
sent to the printer. The output plotter file, IFDOUT 


PLOTTER, serves as an input file for the plotting program. 


D. RUNNING THE PLOTTING PROGRAM 

A Tektronix-618 terminal is used to run the plotting 
program. The first step is to log onto the terminal in the 
normal manner and then define additional storage space by, 


(1) Enter: DEF STOR 1M 
(2) Enter: I CMS 


The plotting program has filename and filetype PLOTIFD 

FORTRAN. To compile the program, 
Enter: FORTGI PLOTIFD 
Unless the program is changed it need only be compiled one 
time. To run the program 
Enters PLOTIFD 

The user will be promptr i for axes and smoothing 
information; enter values and responses as appropriate. The 
transmission loss curve will be displayed on the CRT screen, 
a hard copy may be obtained by pushing the HARD COPY button 
under the screen, and the screen may be cleared by pushing 


the ENTER key. 
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APPENDIX C: PLOTIFD PROGRAM LISTING 


0999000 900000000099090900009300009009909099909990999090 
SYAPH IL DRDHVAY LN OL IAS AN FOOT AAMaMAMAae eee Tee 
9Q99900 00 Catt ted t ttt SANNA NNN OM MM MAM FP HPPPt oe 
+ 3O'JNQD092I9QOU99V99990009000000380090000090009099909099V709 
~“ F QN9990009093999099909099099098099999996999999999999 
=. Q00000000399000000000009900909003900990003900099000S50 
“ Ja AI odd 4S 43 444 JIS I ts J dt Ss J I J 
- GOeaeadaaaaeaa@bAQGAaGagaAaaadrtaaaagaaagaaagagaaagaaaanaaaga 
oe ' 
a et 
Zr 
" =e Ud 
oy ze - 
a. | 3 2 Q 
f Su w us 
a es J 
. be towed 
4,4 Ue 2 a. 
6 [=) < =z 
ue 9 
one 
x a2 
wo z be 
oe Ta =m a 
: re 3 3 ” 
a ‘ 
ae eCQww x 4 
‘i O<e- aad uw 
vs < fod 
= eA ey ud 
WV AQ éK x 
wae 8 rz 
tg VA 
a . 4 & = w 
= 2TtKH< U, b 4 Ww «a 
‘ YtS2 a) f= a 
“Zo baad Qo uy 
Muga e am x «a z 
AGW " uw < (en ] qQ -— ( 
- owe ~ uw = - = 
aa O a (4 hag z z [oa] a 
az A x o—¢ nd j=) Q 
ZO os £ mu o S) 
4 OW "On ~~ «6 i) 2 ” i) 2 fo] =z 
y . rH OL te a 4 — om fered ) a] “4 w“ 
1 w“ a & - —~ VW? bo Oo ie.) 
. Ais ae Ee wa = o™~ CUO ”n < Ww <= Q Q 
: grassy 80g SSZ5Y 9-2 2 Ze 5 Of 
>) on 
A WSSe “e295 © 3 QO 12x 3 a> E roo Q 
SEMQWMWNGaO a] <x<< 3 es 2) tx» oO Us 
4 LOW, x=KH ~ © exzox Oa od > oO vi) 
i, Ce aguonco.a a ~ fw we > - x= = =) 
: rou - Ox Ow’ a La YW - 
i aoe ee es) 3—900m 2 1 Om z» “nz 
y Ndh-~L— B Q - ce = KeEHO XM » ~ & oo aie! 
o a on ra Ai re at fe 2 hon > Q 
4 Omwreo a Oo =e & OKT Zaetiem > Lew Zand= Gg —=o. 
2 AQE co] Om | MIONWOXM ete OND WM Bow Quz 
a <20 Us ° <x wiz ° one 0 a0 
> tee ST ~a * SLI e Ir) a A ort Dares HOW) an- 
‘| Pe veye Ss «< SON URE OF Om PU OW rR OWL ME HL Oe Oe 
! qaw w Us QA N OC ef ETZUO oO OF OO ORONO MNO OF N ZO 
% [odor aay om od xr wy QD WU Swain usth | fw | Zonda Twi 
mat, hote oie b -aa  e Y4- Y De ZO auw~ 
Sy) w Oo ON 4 RO RO FOe RUUROw mr RO 
Sy a <~ CO oO or Tn loo Mi oot tla ro) zm 
a ian = WA 2 ariel he MUCWWCWMNWEWNWCWORWLOAS Occu 
ia an a4 Sd MO RNLGBLwWOSLwMOELORIMSAN CBee 
vw ~ au <« 
ms = =* ~ : * * a % * * * 
hf =~ < # * * + + % 
se - 90 A * *# * + * * + % % 
* ot 
“a ° as 
r VOVVOVVVULWYOUW Y OO VOY Oo Yo ro) re) 7) ro) Qu 


a: 
Wg 103 
§ 

4 


Pe ere eee ee ee ee 
. ". - 


a . a . . =e a - “. . . 
To. matt ee tee ee 


0990990 3900009090999999906999699939009999099000009059 
PDA F i Ole WADA AHO Loopy babe Papi tal cg pad dpe anh aed de A pang he AR 
FRAADAANRIA DOG 000.0 DOOR RMN Misee aD BDODDDD VOKNAHDARDAH 
9990090 90900090909093900090079099000900009909099902N9990199U 
3990999 290029900999960990909099099939393909390909003999909839933A0 
39Q090000990090090000990900000090000000390009990900303900909 
cc mle med me me mS md dd dS me md ed a med nd md dt ad ed med med ad hd) dd 2) nh Sd wed dd dd hd dd) od od 
Che a ho he ha tae Che he, er Ce, Gro hn hy Dh he Chap Te Cg, We hes Ch er Chas he Char Cs Shan Dg Ce (an an Cha, hes he, C- as e, ee hr he, he he Le Ds Ay, he 2, 


WY 
ag 
=) 
Q—4 
uj 
> 
$2 
Sew 
Ww 
J! 
gu 
Woe ud 
Lae] 
J 
Usk 
a2=> ~ 
ar ~ 
<t de ~~ 
eS ~~ ~ 
oO 2 fag N 
1 Gr ~ o 
¢ 
4545 3 age 
OU ie.) On 
Mm 06 Q i 
z =H a | ad 
Uw AQ < eD 
~~ me)p— . oe = 
am ee OWN << e fo] ae 
< ta 2S Wit e ~e 
= bw 1 0 tome pm —4 Sey) fo] or) 
< 27 OC m=zz FD ~ aw 
- a ++ aa, OMmn OZ wa 
rs} w SoSS. vaag a: ouan 
wus wey 
°o So =f a bre actw ~ ONS» 
a ~ + Ga Sos ~ Pee 
+ ~ seas N 
“oO on” aad and wuSin ~ 
a - - = eth, ele Bote ao eeae 
da bored z2iz aaa a RL) ~e ed en | 
wod = Q 3 mW AS —=2 oo ~ Len 
~—~FrOo 3 o + Ore Ne B® wo ore shew & 
ua ~ =z aaexé 2a WWI 9 <g SO ag e we (am 
=a un - Ye == De iB o- 2 00 QxXO 
2323 qa = fee i Wu) aAOWEO wad ON & USL & & 
Dis omg Kk So zt he <ZOer <~ at em Oa, Taya 
me Oct e« QO WOO, O. a > % au eif\ ZOr x & 
met e+ F ma a4 MESES oo O *tONwredWe «oh 
Sw te Ww a | us ome SS SS SS 
dD ad wd w © (6 bbb Up rrm ue Cjue Twwwwwe o Ae 
OQ e« aN Loe ya ri lamane™ pees << AN OD OL 4 OF OD he 9) wr ra BE ews See 
AVA JI zy Zr e ol » HAOMON EMME 
we wo = ame e Owns =D matt tl VCWMNZOOULTLOW > 
Wwe Wy WY VEO Wa@mr WwW i] CS WAIL ws ren org peg ak ae KY 
<< oa “aH a | ome USE EO ibe 2 tua 
UW, Oo 22 20 oOo < mat ZR OO Ee Me we TOU 
[a fe a) iS aI Oo WY Qrwwrdt 
ld [volo We 4 eee ee ee ee ee eee ee ee 
* =z * ReRRH Pe re oe ee ee ee el 
- Ou + seit Q QO aA44¢ oda q 
Um aed Q VQVOQUQQOUVIVOOUO 


0990090000000000909009 
RAROANY LAS VDRO PY 
ARKOIO QIOQ00N~e 
9] QEDO ado cred mad 0k pend a a eh rk poe ad Ob rad ad red 
% : 990009909099999009009 
ae * 99000000009000e8990009 
4 cel cd aed ad ah med ed ad ded aed ed ad ad nd) edd = 
4 0.2.3.3.40332.0.36244638.234940. 


@ 
8) 
IN 
Y 


MISSIGN LOSS CURVE RANGES*,/ 
3 OR coo 1°) 


Noe 


x 
A 
D e 
The “Olt 
YOU WANT AVERAGE: aT EACH STeP TO *, 
( ENTERS 2 OR 


e 
tOw MANY POINTS DO 


MCOTHING??,/ 5° 


a, “@q&Qee fee we 


2u ee dl ee 

[o4 = = 

aad Setisssrossro 

SS & @- +" Se-- Seo" s 
eure + * 


° ANIM THO 
° aoo 3900 
dl a BAD ABH 


105 


a NE an ne oe on oi Se Se SS 
ap a To Ww Fs he ea ee srs 


APPENDIX D 
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